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Abstract 

The  ever-increasing  demand  in  surveillance  is  to  produce  highly  accurate  target  and  track  identifica¬ 
tion  and  estimation  in  real-time,  even  for  dense  target  scenarios  and  in  regions  of  high  track  contention. 
The  use  of  multiple  sensors,  through  more  varied  information,  has  the  potential  to  greatly  enhance  target 
identification  and  state  estimation.  For  multitarget  tracking,  the  processing  of  multiple  scans  all  at  once 
yields  the  desired  track  identification  and  accurate  state  estimation;  however,  one  must  solve  an  NP-hard 
data  association  problem  of  partitioning  observations  into  tracks  and  false  alarms  in  real-time.  This  report 
summarizes  the  development  of  a  multisensor-multitarget  tracker  based  on  the  use  of  near-optimal  and  real¬ 
time  algorithms  for  the  data  association  problem  and  is  divided  into  several  parts.  The  first  part  addresses 
the  formulation  of  multisensor  and  multiscan  processing  of  the  data  association  problem  as  a  combinatorial 
optimization  problem.  The  new  algorithms  under  development  for  this  NP-hard  problem  are  based  on  a 
recursive  Lagrangian  relaxation  scheme,  construct  near-optimal  solutions  in  real-time,  and  use  a  variety  of 
techniques  such  as  two  dimensional  assignment  algorithms,  a  bundle  trust  region  method  for  the  nonsmooth 
optimization,  graph  theoretic  algorithms  for  problem  decomposition,  and  a  branch  and  bound  technique  for 
small  solution  components.  A  brief  computational  complexity  analysis  as  well  as  a  comparison  with  some 
additional  heuristic  and  optimal  algorithms  is  included  to  demonstrate  the  efficiency  of  the  algorithms.  A 
two  radar  system  with  data  supplied  by  Rome  Labs  is  used  to  demonstrate  the  efficiency  and  robustness  of 
current  multisensor-multitarget  tracker  that  is  based  on  these  fast  data  association  algorithms. 
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1.  Introduction.  The  ever-increasing  demand  in  surveillance  is  to  produce  highly  accurate  target  and 
track  identification  and  estimation  in  real-time,  even  for  dense  target  scenarios  and  in  regions  of  high  track 
contention.  The  use  of  multiple  sensors,  through  more  varied  information,  has  the  potential  to  greatly  improve 
state  estimation  and  track  identification.  This  approach  is  part  of  a  much  broader  problem  called  data  fusion, 
which  for  military  applications  is  defined  as  “a  multilevel,  multifaceted  process  dealing  with  the  detection, 
association,  correlation,  estimation  and  combination  of  data  and  information  from  multiple  sources  to  achieve 
refined  state  and  identify  estimation,  and  complete  and  timely  assessments  of  situation  and  threat”  [64].  The 
various  problems  are  generally  partitioned  into  three  or  more  levels:  (1)  fused  position  (state)  and  identity,  (2) 
hostile  or  friendly  military  situation  assessments,  and  (3)  hostile  force  threat  assessments.  (Comprehensive 
discussions  can  be  found  in  the  books  of  Waltz  and  Llinas  [64]  and  Hall  [26].)  Level  1  deals  with  single  and 
multisource  information  involving  tracking,  correlation,  alignment,  and  association  by  sampling  the  external 
environment  with  multiple  sensors  and  exploiting  other  available  sources.  Numerical  processes  thus  dominate 
Level  1;  symbolic  reasoning  involving  various  techniques  from  artificial  intelligence  permeate  Levels  2  and 
3.  This  report  focuses  on  Level  1  data  fusion  with  the  goal  being  to  use  multiple  sensors  to  achieve  superior 
state  estimation  and  track  identification. 

Sensor  fusion  systems  vary  greatly  depending  on  the  particular  needs  of  a  surveillance  system.  Key 
issues  in  the  design  of  such  a  system  include  sensor  type  (active  or  passive),  sensor  location  (distributed 
or  collocated),  and  the  level  of  association,  which  ranges  from  sensor  level  fusion  to  centralized  fusion  with 
hybrids  in  between.  Although  there  are  many  such  issues,  the  central  problem  in  any  surveillance  system  is 
the  data  association  problem  of  partitioning  measurements  into  tracks  and  false  alarms.  To  explain  this  data 
association  problem,  we  must  first  address  the  levels  of  association. 

At  one  extreme  is  sensor  level  tracking,  wherein  each  sensor  forms  tracks  from  its  own  measurements  and 
then  the  tracks  from  the  sensors  are  fused  in  a  central  location.  Once  the  correlation  is  complete,  one  then 
combines  the  tracks  with  appropriate  modification  in  the  statistics  [13].  Compared  with  central  level  fusion, 
the  advantages  include  the  reduced  communication  costs  between  the  sensors  and  central  processing  unit 
and  easier  data  association.  The  disadvantages  are  that  combined  track  estimates  tend  to  be  worse  than  in 
central-level  fusion  and  the  error  independence  assumptions  in  data  association  are  no  longer  valid,  thereby 
introducing  additional  complexity  into  the  problem  [11,13].  At  the  other  extreme  is  centralized  fusion  in 
which  sensors  send  measurements  to  a  central  processing  unit  where  they  can  be  combined  to  give  superior 
state  estimation  [11]  (compared  to  fusion  of  sensor  level  tracks).  The  difficulties  are  generally  claimed  to  be 
data  association  ( our  strength ),  communication  costs  between  the  sensor  and  central  processing  unit,  and 
the  loss  of  the  tracking  capability  if  the  central  processing  unit  becomes  inoperative.  In  reality,  current 
and  proposed  sensor  fusion  systems  for  any  surveillance  system  make  use  of  both  systems.  Certainly,  one 
can  treat  a  hybrid  of  these  two  systems  by  sending  the  observations  associated  with  a  track  obtained  at 
the  sensor  level  to  a  central  processing  unit  and  treat  the  association  as  in  centralized  fusion  [11].  Having 
explained  the  level  of  data  association,  we  now  return  to  a  brief  overview  of  the  methods  of  data  association 
for  central  and  some  hybrid  central-sensor  level  tracking. 

The  existing  algorithms  range  from  single  scan  or  sequential  processing  to  multiscan  processing.  Methods 
for  the  former  include  nearest  neighbor,  one-to-one  or  few-to-one  assignments,  and  all-to-one  assignments  as 
in  the  joint  probabilistic  data  association  (JPDA)  [6]  in  single  sensor  tracking.  Problems  involving  one-to- 
one  or  few-to-one  assignments  are  generally  formulated  as  (two  dimensional)  assignment  or  multi-assignment 
problems  for  which  there  are  some  excellent  algorithms  [7,  8,  9].  This  methodology  is  real-time  but  can  result 
in  a  large  number  of  partial  and  incorrect  assignments,  particularly  in  dense  or  high  contention  scenarios, 
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and  thus  incorrect  track  identification.  The  difficulty  is  that  decisions,  once  made,  are  irrevocable,  so  that 
there  is  no  mechanism  to  correct  misassociations.  The  use  of  all  observations  in  a  scan  (e.g.,  JPDA)  to 
update  a  track  moderates  the  misassociation  problem  and  has  been  successful  for  tracking  a  few  targets  in 
dense  clutter  [6]. 

Deferred  logic  techniques  consider  several  data  sets  or  scans  of  data  from  multiple  sensors  all  at  once 
in  making  data  association  decisions.  At  one  extreme  is  batch  processing  in  which  all  observations  (from 
all  time)  are  processed  together,  but  this  is  computationally  too  intensive  for  real-time  applications.  The 
other  extreme  is  sequential  processing.  Deferred  logic  methods  between  these  two  extremes  are  of  primary 
interest  in  this  work.  The  key  advantage  of  this  approach  is  the  ability  to  change  data  association  decisions 
over  several  of  the  most  recent  scans  of  data.  It  is  this  feature  that  leads  to  superior  track  estimation .  The 
principal  deferred  logic  method  used  to  track  large  numbers  of  targets  in  low  to  moderate  clutter  is  called 
multiple  hypothesis  tracking  (MHT)  in  which  one  builds  a  tree  of  possibilities,  assigns  a  likelihood  score  based 
on  Bayesian  estimation,  develops  an  intricate  pruning  logic,  and  then  solves  the  data  association  problem  by 
explicit  enumeration  schemes.  The  fundamental  limitation  of  MHT,  as  it  now  exists,  is  that  it  is  an  NP-hard 
combinatorial  optimization  problem,  so  that  in  dense  scenarios  and  high  track  contention  or  with  multiple 
sensor  input,  the  time  required  to  solve  this  problem  optimally  can  grow  exponentially  with  the  size  of  the 
problem.  This  failure  is  not  graceful,  i.e.,  the  method  is  not  robust  with  respect  to  real-time  needs.  Thus 
to  make  MHT  viable,  near-optimal  algorithms  are  needed  to  solve  the  data  association  problems  to  the  noise 
level  in  real-time.  This  is  precisely  the  subject  of  this  research  program  and  report.  The  centralized  fusion 
approach  is  also  highly  parallelizable  and  is  ripe  for  the  use  of  parallel  computers  in  the  future. 

The  first  topic  (Section  2)  in  this  report  is  the  formulation  of  multisensor  and  multiscan  processing  of  the 
data  association  problem  as  an  NP-hard  combinatorial  optimization  problem.  Next,  an  overview  of  some  of 
the  near-optimal  and  real-time  algorithms  for  solving  this  problem  is  presented  in  Section  3.  The  algorithms 
under  development  are  based  on  a  recursive  Lagrangian  relaxation  scheme,  construct  near-optimal  solutions 
in  real-time,  and  use  a  variety  of  techniques  ranging  from  two  dimensional  assignment  algorithms,  a  bundle 
trust  region  method  for  the  nonsmooth  optimization,  graph  theoretic  properties  for  problem  decomposition, 
and  a  branch  and  bound  technique  for  small  solution  components.  These  topics  are  presented  in  Sections  3 
and  4.  A  brief  computational  complexity  analysis  as  well  as  a  comparison  with  some  additional  heuristic  and 
optimal  algorithms  is  included  in  Section  5  to  demonstrate  the  efficiency  of  the  algorithms.  The  existing  and 
on-going  software  work  is  discussed  in  Section  6.  Finally,  in  Section  7,  a  two  radar  system  with  data  supplied 
by  Rome  Labs  is  used  to  demonstrate  the  efficiency  and  robustness  of  current  multisensor-multitarget  tracker 
that  is  based  on  these  fast  data  association  algorithms. 

2.  Formulation  of  the  Data  Association  Problem.  The  goal  of  this  section  is  to  explain  the  formulation 
of  the  data  association  problem  that  governs  large  classes  of  data  association  problems  in  centralized  or 
hybrid  centralized-sensor  level  multisensor /multitarget  tracking.  The  presentation  is  brief;  technical  details 
are  presented  for  both  track  initiation  and  maintenance  in  [49]  for  nonmaneuvering  targets  and  [52]  for 
maneuvering  targets.  These  works  also  contain  expressions  for  the  likelihood  ratios  Li1...iN  used  in  the 
score  in  the  following  equations  (2.3)  and  (2.4).  The  formulation  presented  here  is  of  sufficient  generality 
to  cover  the  MHT  work  of  Reid  [61],  Blackman  and  Stein  [10],  and  modifications  by  Kurien  [34]  to  include 
maneuvering  targets.  As  suggested  by  Blackman  [10]  this  formulation  can  also  be  modified  to  include  target 
features  (e.g.,  size  and  type)  into  the  scoring  function. 

The  data  association  problems  for  multisensor  and  multitarget  tracking  considered  in  this  work  are 
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generally  posed  [6,10,35,49,52]  as  that  of  maximizing  the  posterior  probability  of  the  surveillance  region 
(given  the  data)  according  to 


Maximize 


p(r  =  7  i  zN) 

P(T  =  7°  |  ZN ) 


(2.1) 


where  ZN  represents  N  data  sets,  7  is  a  partition  of  indices  of  the  data  (and  thus  induces  a  partition  of  the 
data),  T*  is  the  finite  collection  of  all  such  partitions,  T  is  a  discrete  random  element  defined  on  T*,  70  is 
a  reference  partition,  and  P(T  —  7  |  ZN )  is  the  posterior  probability  of  a  partition  7  being  true  given  the 
data  ZN .  The  term  partition  is  defined  below;  however,  this  framework  is  currently  sufficiently  general  to 
cover  set  packings  and  coverings  [35]. 

Consider  N  data  sets  Z(k)  (k  =  1, . . . ,  N)  each  of  Mk  reports  {zfk  ,  and  let  ZN  denote  the  cumulative 
data  set  defined  by 

m  =  {4X=i  and  Zn  =  {Z(1),...,Z(N)},  (2.2) 


respectively.  In  multisensor  data  fusion  and  multitarget  tracking  the  data  sets  Z(k)  may  represent  different 
classes  of  objects,  and  each  data  set  can  arise  from  different  sensors.  For  track  initiation  the  objects  are 
measurements  that  must  be  partitioned  into  tracks  and  false  alarms.  In  our  formulation  of  track  maintenance 
[49,52],  which  uses  a  moving  window  over  time,  one  data  set  will  be  tracks  and  remaining  data  sets  will  be 
measurements  which  are  assigned  to  existing  tracks,  as  false  measurements,  or  are  assigned  to  initiating 
tracks.  In  sensor  level  tracking,  the  objects  to  be  fused  are  tracks  [10].  In  centralized  fusion  [10],  the  objects 
may  all  be  measurements  that  represent  targets  or  false  reports,  and  the  problem  is  to  determine  which 
measurements  emanate  from  a  common  source. 

We  specialize  the  problem  to  the  case  of  set  partitioning  [49]  defined  in  the  following  way.  First,  for 
notational  convenience  in  representing  tracks,  we  add  a  dummy  report  z *  to  each  of  the  data  sets  Z(k)  in 
(2.2),  and  define  a  “track  of  data”  as  (z/  , . . . ,  z(*N)  where  4  and  z\k  can  now  assume  the  values  of  0  and  zfi, 
respectively.  A  partition  of  the  data  will  refer  to  a  collection  of  tracks  of  data  wherein  each  report  occurs 
exactly  once  in  one  of  the  tracks  of  data  and  such  that  all  data  is  used  up;  the  occurrence  of  a  dummy  report 
is  unrestricted.  The  dummy  report  Zq  serves  several  purposes  in  the  representation  of  missing  data,  false 
reports,  initiating  tracks,  and  terminating  tracks  [49].  The  reference  partition  is  that  in  which  all  reports 
are  declared  to  be  false. 

Next  under  appropriate  independence  assumptions  one  can  show  [49] 


p(r  =  7  |  zN)  _  _ 

P(T  =  70  |  ZN)  ~  7  ~ 


(2.3) 


L{l...iN  is  a  likelihood  ratio  containing  probabilities  for  detection,  maneuvers,  and  termination  as  well  as 
probability  density  functions  for  measurement  errors,  track  initiation  and  termination.  Then  if  = 

In  L{  1  ■  ■  ■  t  jv  1 


—  In 


-piMz^y 

.P(l°\ZN). 


=  E 


Ci 


*l*~* N * 


*tv)67 


(2.4) 


Using  (2.3)  and  the  zero-one  variable  Zil...{N  =  1  if  (?, , . .  .,ijv)  6  7  and  0  otherwise,  one  can  then  write  the 


5 


problem  (2.1)  as  the  following  N-dimensional  assignment  problem: 

Mi  M  jv 

Minimize  Y2  '  ^2 

t  i =0  tV=o 

M2  M  jv 

Subject  To:  ^  ^  =  1,  h  = 

*2  =  0  *AT=0 

Mi  Mk-i  Mk~ |-i  Mat 

V  V  V  (2.5) 

*1=0  U-I=0*fc+1=0  *v  =  0 

for  4  =  1,...,  Affc  and  k  =  2, . . . ,  AT  —  1, 

Mi  Mjv-i 

E  ■  •  ■  E  =  ^  ijv  =  1,  ...,Afjv 

*1=0  *V-i  =0 

z»! •••«„  G  {0,  !}  f°r  a11  U,  •  •  • .  *jv, 

where  co...o  is  arbitrarily  defined  to  be  zero.  Here,  each  group  of  sums  in  the  constraints  represents  the  fact 
that  each  non-dummy  report  occurs  exactly  once  in  a  “track  of  data” .  One  can  modify  this  formulation  to 
include  multiassignments  of  one,  some,  or  all  the  actual  reports.  The  assignment  problem  (2.5)  is  changed 
accordingly.  For  example,  if  z\k  is  to  be  assigned  no  more  than,  exactly,  or  no  less  than  n\k  times,  then  the 
“  =  1”  in  the  constraint  (2.5)  is  changed  to  “  <,  =,  >  n\k ,  ”  respectively.  Modifications  for  group  tracking 
and  multiresolution  features  of  the  surveillance  region  will  be  addressed  in  future  work.  In  making  these 
changes,  one  must  pay  careful  attention  to  the  independence  assumptions  that  need  not  be  valid  in  many 
applications. 

For  track  maintenance ,  we  use  a  sliding  window  of  N  data  sets  and  one  data  set  containing  established 
tracks  [49,52].  The  formulation  is  the  same  as  above  except  that  the  dimension  of  the  assignment  problem 
is  now  N  +  1. 

3.  Overview  of  the  Lagrangian  relaxation  algorithms.  Having  formulated  an  N-dimensional  assign¬ 
ment  problem  (2.5),  we  now  turn  to  a  description  of  the  Lagrangian  relaxation  algorithms.  The  algorithms 
presented  here  are  the  subject  of  numerous  publications  [42,  45,  46,  47,  50,  55,  59]  and  two  patents  [53,58]  and 
were  developed  under  the  current  and  previous  AFOSR  grants.  The  relaxation  procedures  presented  here  are 
based  on  relaxing  an  arbitrary  number  of  constraints.  Thus,  subsection  3.1  presents  many  of  the  relaxation 
properties  associated  with  the  relaxation  of  an  n-dimensional  assignment  problem  to  an  m-dimensional  one 
via  a  Lagrangian  relaxation  of  n  —  m  sets  of  constraints.  Although  any  n  —  m  sets  can  be  relaxed  the  descrip¬ 
tion  here  is  based  on  relaxing  the  last  n  —  m  sets  of  constraints  and  keeping  the  first  m  sets.  Given  either  an 
optimal  or  suboptimal  solution  of  the  relaxed  problem,  a  technique  for  restoring  feasibility  to  the  original  n- 
dimensional  problem  is  presented  in  subsection  3.2  and  an  overview  of  the  Lagrangian  relaxation  algorithm, 
in  subsection  3.3.  A  presentation  of  these  various  algorithms  is  the  subject  of  forthcoming  publications  [59, 
59,  55]  and  the  thesis  of  Alexander  J.  Robertson  [60]. 

The  following  notation  will  be  used  throughout  the  remainder  of  the  work.  Let  N  be  an  integer  such 
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that  N  >  3  and  let  n£{3,.. . ,  N}.  The  n-dimensional  assignment  problem  is 

Mn 

Minimize  Vn(z)  =  £  ' ' '  £  tfi-fX-i. 

*1=0  Irt=0 


m2 

Mn 

E 

•••E  <•••*.= 

1,  ii  =  1, . . . ,  Mi, 

i  2  =0 

*  n  — 0 

Mi 

Mk-i  Mfe+1 

Mn 

E 

-  E  E 

*1=0 

*fe  —  i =0 tfc  +  i =o 

*  n  —  0 

for  4  =  1,..., 

Mk  and  k  =  2, . . . ,  n  — 

Mi 

M„-i 

E 

...  T  zn  . 

II 

i—i 

. 

3 

II 

£ 

*1=0 

*  n  —  1  =  0 

(3.1) 


E  {0, 1}  for  all  ij, . . . ,  in. 


To  ensure  that  a  feasible  solution  of  (3.1)  always  exists,  all  variables  with  exactly  one  nonzero  index  (i.e., 
variables  of  the  form  ^o—o*fco— o  for  ^  are  assumed  free  to  be  assigned  and  the  corresponding  cost 
coefficients  are  well-defined. 

3.1  The  Lagrangian  Relaxed  Assignment  Problem.  The  n-dimensional  assignment  problem  (3.1)  has 
n  sets  of  constraints.  A  (M*  +  l)-dimensional  multiplier  vector  associated  with  the  k- th  constraint  set  will  be 
denoted  by  uk  =  («§ ,  it* , . . . ,  uhMk)T  with  uk  =  0  and  k  —  1, . . . ,  n.  The  n-dimensional  assignment  problem 
(3.1)  is  relaxed  to  an  m-dimensional  assignment  problem  by  incorporating  n  —  m  of  the  n  sets  of  constraints 
into  the  objective  function.  Although  any  constraint  set  can  be  relaxed,  the  description  of  the  relaxation 
procedure  for  (3.1)  will  be  based  on  the  relaxation  of  the  last  n  —  m  sets  of  constraints.  The  relaxed  problem 
is 

Mi  Mn 

Minimize  ?W(zn;  um+1, . . . ,  un)  =  ^’ ' '  £  c"i •  •••«» 

*1=0  t'rt=0 

n  Mk  r  Mk~i  Mk+ 1  Mn 

+  E  E  <  E  E  E  •••EC-..-1 

*i—0  i  =0  * =0  * n  =0 

Mi  M  n  n  n  Mk 

sE  -E [«?,..+  E  E  E< 

*1=0  *  n — 0  k—m+ 1  k=m+ 1  ifc= 0 

M2  Mn 

Subj.To  X)  •  •  •  =  1,  *1  = 

*2  =  0  *n  =  0 

Mi  Mk- 1  Mk+i  Mn 

£•••£  £  •••£-%•  ,„  =  i- 

*1=0  *fe_i  =0  ik+i=0  * n=0 

for  ik  =  1, . . . ,  Mk  and  k  =  2, . . . ,  m, 

€  {0, 1}  for  all  n, . . . ,  in- 

An  optimal  (or  suboptimal)  solution  of  (3.2)  can  be  constructed  from  that  of  an  m-dimensional  assign¬ 
ment  problem.  To  show  this,  define  for  each  (ii,...,im)  an  index  (jm+ii  •  •  • ,  jn)  — 


(3.2) 
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■  ■  ■  ,im),  ■  ■  ■  ,jn(iu  ■  ■  ■  ,im))  and  a  new  cost  function  by 


(im+i,  ■  •  -  ,jn)  =  arg  min<  ^  u. 


k—m  +  l 
n 


4  =  0,...,  Mk  and  k  =  m  +  1, . . . ,  n 


c™-im  =  c"i +  2  uf* for (*i» •  •  •  ■  *■») # (°> •  •  •  > o) 


(3.3) 


Mm  +  i  Af  n  ✓  n  s 

=  E  ---E  Minimum]  0>cg...Wm+1...im+  E  A  } 

im+1=0  in—Q  ^  k=m+ 1  ' 


Cq-o 


(If  (jm+i , . . . ,  jn)  is  not  unique,  choose  the  smallest  such  jm+ 1 ,  amongst  those  (n  —  ra)-tuples  with  the  same 
jm+ 1  choose  the  smallest  jm+2>  etc.,  so  that  (jm+ i,  *  •  .,in)  is  uniquely  defined.)  Using  the  cost  coefficients 
defined  in  this  way,  the  following  m-dimensional  assignment  problem  is  obtained: 


Mi 

Mm 

1.  7/m+ 1 
>  u  >  ■  •  •  > 

Un)=Vm(zm)  = 

;£ 

•£ 

>1=0 

—  0 

m2 

Mm 

£- 

£  C-  i. 

=  1,  *1=1.--- 

Afi, 

1*2  —  0 

l*m  =  0 

Ml 

Mk- 1  Mk+i 

Mm 

£  • 

£  £ 

£  C-.. 

=  1, 

i‘i=0 

ifc_i=0ifc+i  =  i 

o 

II 

•i 

o 

for 

4  =  1, . . . 

. ,  Mk  and  k  =  2, 

. . . ,  m 

-1, 

Mi 

Mm-  1 

£  • 

£  C-, 

im  =  ^  ^ 

. . . ,  Mm 

1*1=0 

zm  . 
*!•••*»» 

6  {0, 1}  for 

all  4 , . . . ,  im . 

(3.4) 


As  an  aside,  observe  that  any  feasible  solution  zn  of  (3.1)  yields  a  feasible  solution  zm  of  (3.4)  via  the 
construction 

zm  m  _  /  1  if  =  1  for  some  (im+ 1, . . . ,  tn) 

11  "*m  lo  otherwise. 

Thus  the  m-dimensional  assignment  problem  (3.4)  has  at  least  as  many  feasible  solutions  of  the  constraints 
as  the  original  problem  (3.1). 

The  following  theorem  states  that  an  optimal  solution  of  (3.2)  can  be  computed  from  that  of  (3.4). 
The  converse  is  contained  in  Theorem  3.2.  Furthermore,  if  the  solution  of  either  of  these  two  problems  is 
e-optimal,  then  so  is  the  other. 

Theorem  3.1.  Let  wm  be  a  feasible  solution  to  problem  (3.4)  and  define  wn  by 


<• 

=  if  (*m  +  l,  •  •  • 

4»)  —  0*m+l  3  •  ■  • 

>  in) 

and  (4, . . 

•  j  1  m  )  7“ 

(o,.. 

-.0) 

<• 

■■>•„  =  0 

if  (^m+l  >  •  • 

•  *n)  7^  C?m+1  j  •  • 

•  >  in  ) 

and  (4, . . 

•  )  z  m )  7^ 

(0... 

...0) 

<-o,m+1. 

=  1 

if  c£-.o,m+I. 

■■«.+  t  < 

<  o 

k=m  +  l 

••>•„  =  0 

lf  c0 -0tm+1- 

E  < 

>  0 

k=m+l 
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Then  wn  is  a  feasible  solution  of  the  Lagrangian  relaxed  problem  (3.2)  and  um+l, . .  - ,  un)  = 

um+l , . . . ,  un )  -  SlUm+i  2of=o  uik  *  ^  2n  addition,  wm  is  optimal  for  (3.4),  then  wn  is  an  optimal 
solution  of  (3.2)  and  $m„(um+1  =  $mr>(um+1  ,...,un)~  ££=m+1  T.it=o  ui„ ■ 


With  the  exception  of  one  equality  being  converted  to  an  inequality,  the  following  theorem  is  a  converse 
of  this  theorem. 

Theorem  3.2.  Let  wn  be  a  feasible  solution  to  problem  (3.2)  and  define  wm  by 
Mm  + 1  Mn 

wh-im=  •••£  K-in  for  (*1.  •■•>*'"»)  #  (°»  •  -  »  °) 

+  l  =  0  i  n~ 0 

n 

w£..0  =  0  if  =  (0.....0)  and  c$...oim+i...in  +  ^  «,•*  >  0  for  all  (iro+i,  •  ■  • ,  *n)  (3-6) 

k=m-\-l 

n 

Wo'-o  =  1  if  (*i>  •••,**»)  =  (0.....0)  and  cg...0,m+1...jn  +  ^  ukik  <  0  for  some  (im+i, . . . ,  »„). 

fe=m+ 1 

Then  wm  is  a  feasible  solution  of  the  problem  (3.4)  and  c/)mn(wn ;  um+1 , . . . ,  un)  >  <f>mn(wm ;  um+l , . . . ,  un )  — 
Y^k=m+i  uik‘  Jn  addition,  wn  is  optimal  for  (3.2),  then  wm  is  an  optimal  solution  of  (3.4), 

. •")  =  #»»(• »";o”+I . «“)  -  E»-m+iE£io<.  . «")  = 

. «")-EL„+,E£o«t 


3.2  The  recovery  procedure.  The  next  objective  is  to  explain  a  recovery  procedure,  i.e.,  given  a  feasible 
(optimal  or  suboptimal)  solution  wm  of  (3.4)  (or  wn  of  (3.2)  constructed  via  Theorem  3.1),  generate  a  feasible 
solution  zn  of  (3.1)  which  is  close  to  wm  in  a  sense  to  be  specified.  We  first  assume  that  no  variables  in  (3.1) 
are  preassigned  to  zero;  this  assumption  will  be  removed  shortly.  The  difficulty  with  the  solution  wn  is  that 
it  need  not  satisfy  the  last  n  —  m  sets  of  constraints  in  (3.1).  (Note  however  that  if  wm  is  an  optimal  solution 
for  (3.4)  and  wn  (constructed  as  in  Theorem  3.1)  satisfies  the  relaxed  constraints,  then  wn  is  optimal  for 
(3.1).)  The  recovery  procedure  described  here  is  designed  to  preserve  the  0-1  character  of  the  solution  wm 
of  (3.4)  as  far  as  possible:  If  =  1  and  ii  ^  0  for  at  least  one  l  =  !,...,  m,  the  corresponding  feasible 

solution  zn  of  (3.1)  is  constructed  so  that  z"  ,-n  =  1  for  some  (im+i,  •  • .  ,in)*  By  this  reasoning,  variables 
of  the  form  zfi '...o*m+l...in  can  ^signed  a  value  of  one  in  the  recovery  problem  only  if  ,0  =  1.  However, 
variables  Zo---oim+i---*n  wdl  ke  treated  differently  in  the  recovery  procedure  in  that  they  can  be  assigned  0 
or  1  independent  of  the  value  WqTq.  This  increases  the  feasible  set  of  the  recovery  problem,  leading  to  a 
potentially  better  solution. 


Let  {(z2, . . .  ,iJm) }j=i  be  an  enumeration  of  indices  of  wm  (or  the  first  m  indices  of  wn  constructed  in 
Theorem  3.1)  such  that  w™  .j  —  1  and  (z{, . . . ,  z£j  ^  (0, . . . ,  0).  Set  (zj, . . . ,  z'^)  =  (0, . . . ,  0)  for  j  =  0  and 
define 


n-m  +  l 

jim+l  -iv. 


C.i 


1 ' '  ’*  m*  m  +  l  "  '*  n 


for  4  =  0,...,  Mk\  k  =  m  +  1, . .  .  ,n;  j  =  0, . . . ,  M0. 


(3.7) 
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Let  Y  denote  the  solution  of  the  (n  —  m  +  l)-dimensional  assignment  problem 

Mo  + 1  Mn 

Minimize  E  E  ’ "  E 

j=0  im  +  i=0  *n  —0 

Mm+i  Mn 

Subject  To  E  '  • '  E  yjim+i-in  =  1.  j  =  l,...,M0, 

tm+l=0  *tt=0 

M0  Mm+ 1  Mk-i  Mfc+i  Mn 

EE  T:  E  -E  =  1.  <3-«> 

;=0  t'm+i=0  u_i=0  u+i=0  *n=0 

for  4  =  1,...,  Mjfc  and  fc  =  m  +  1, . . . ,  n  —  1, 

Mo  M  m+ 1  Mn_i 

E  =  =  1)  •  •  •  !  -^nj 

i=0»'m+i=0  »n-l=0 

yjim+i-in  e  {0,1}  for  all  ■••>*»• 

The  recovered  feasible  solution  z"  of  (3.1)  corresponding  to  the  multiplier  set  {um, . . . ,  u"}  is  then  defined 
by 

z?  _/l  if  =  (*!.•••,  *L)  for  some  j  =  0,...,  Mo  and  Yjim+l...in  =  1  (39) 

81  ",n  lO  otherwise. 

This  recovery  procedure  is  valid  as  long  as  all  cost  coefficients  cn  are  defined  and  all  zero-one  variables  in  zn 
are  free  to  be  assigned.  Modifications  are  necessary  for  sparse  problems.  If  the  cost  coefficient  c”  .j  . 

*  1  ‘ '  '*  m*  m-fl  '  ‘ -t  « 


is  well  defined  and  the  zero-one  variable  zn2 


n-m  +  l  __  rn 

j*m+l‘‘*n  *j— <» 


as  in  (3.7)  with  z™{ 


®m*m+l  ”’*n 

n-m+l 

m-f-1  ‘ '  ‘*n 


is  free  to  be  assigned  to  zero  or  one,  then  define 
being  free  to  be  assigned.  Otherwise, 


n-m  +  l 


IS 


preassigned  to  zero  or  the  corresponding  arc  is  not  allowed  in  the  feasible  set  of  arcs.  To  ensure  that  a  feasible 
solution  exists,  we  now  only  need  ensure  that  the  variables  are  free  for  j  =  0,1,...,  Mo.  (Recall  that 

all  variables  of  the  form  2o-.-*v*-o  are  free  (f°  be  assigned)  and  all  coefficients  of  the  form  are  well 

defined  for  k  =  1, . . . ,  n.)  This  is  accomplished  as  follows:  If  the  cost  coefficient  e”  .3 


nO  -O 


is  well  defined  and 


znj-j  j  is  free,  then  define  c”0..™+1  =  c”.,  with  being  free.  Otherwise,  since  all  variables 

of  the  form  ztf. ..ifc ...0  are  known  feasible  and  have  well-defined  costs,  put  Cjq.™+1  =  ^jT=i*j^o  co-..o*io  -o4 
3.3.  Summary  of  the  Lagrangian  Relaxation  Algorithm.  Starting  with  the  N- dimensional  assign¬ 
ment  problem  (3.1),  i.e.  n  =  N,  the  algorithm  described  below  is  recursive  in  that  the  TV-dimensional 
assignment  problem  is  relaxed  to  an  m-dimensional  one  by  incorporating  (n  —  m)  sets  of  constraints  into 
the  objective  function  using  a  Lagrangian  relaxation  of  this  set.  This  problem  is  maximized  with  respect 
to  the  Lagrange  multipliers,  and  a  good  suboptimal  solution  to  the  original  problem  is  recovered  using  an 
(n  —  m  +  l)-dimensional  assignment  problem.  Each  of  these  two  (the  m-dimensional  and  the  (n  -  m  +  1)- 
dimensional  assignment  problems)  can  be  solved  in  a  similar  manner.  Here  we  describe  one  loop  in  this 
procedure. 


The  Lagrangian  Relaxation  Algorithm  for  the  n-Dimensional  Assignment  Problem 

Assume  TV  >  3  and  choose  n  €  {3, . . .,  N}  and  2  <  m  <  n.  To  obtain  a  near-optimal  solution  of  the 
n-dimensional  assignment  problem  (3.1),  proceed  as  follows: 

A.  Initialization:  Choose  an  initial  approximation  {u™+1, . . . ,  Uq}. 

B.  Use  a  non-smooth  optimization  technique  (see  subections  4.2  and  4.3)  to  maximize  <£mn  in  (3.2), 
i.e., 

Maximize  {<£mn(itm+1 , . . . ,  un)  |  uk  E  IRMfc+1  ;  k  =  m  +  1, . . . ,  n}  (3.10) 
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C.  Given  an  optimal  solution  (um+1, . . .,  un)  of  (3.10),  compute  an  optimal  solution  of  (3.4). 

D.  Using  the  procedure  described  in  Section  3.2,  solve  the  (n  — m+  l)-dimensional  assignment  problem 
(3.8)  to  recover  a  feasible  solution  zn  of  the  n-dimensional  assignment  problem  (3.1). 

3.4.  Comments  on  the  Various  Algorithms. 

Of  the  many  algorithms  descibed  above  we  dicuss  them  as  three  different  classes  of  algorithms.  The 
first  is  from  our  earlier  work  [42,  50]  and  is  that  in  which  one  set  of  constraints  is  relaxed,  yielding  an 
m  =  n  —  1  dimensional  problem  by  incorporating  one  set  of  constraints  into  the  objective  function  via 
the  Lagrangian.  We  maximize  the  relaxed  problem  with  respect  to  the  corresponding  Lagrange  multipliers 
and  then  reconstruct  a  feasible  solution  to  the  n-dimensional  problem  using  a  two-dimensional  assignment 
problem.  This  first  relaxation  scheme  describes  the  framework  for  our  earlier  algorithm  work  [42,  50]  and 
also  incorporates  the  algorithm  presented  in  [15].  The  algorithm  of  Deb  et.al.  is  a  special  case  of  ours  in 
that  only  one  nonsmooth  iteration  is  taken  at  each  inner  relaxed  m-dimensional  level,  no  merit  function 
is  implemented,  and  no  decomposition  is  performed.  Their  nonsmooth  optimization  method  is,  however, 
different  from  those  that  we  use. 

The  second  algorithm  is  somewhat  of  a  mirror  image  of  the  first  in  that  n  —  2  sets  of  constraints  are 
relaxed,  yielding  an  m  =  2  dimensional  problem.  A  feasible  solution  to  the  n  dimensional  problem  is 
then  recovered  using  an  n  —  1  dimensional  problem.  In  this  case  the  function  values  and  subgradients  of 
$2n(u3, . . .,  un)  can  be  computed  optimally  via  a  two  dimensional  assignment  problem.  The  significant  ad¬ 
vantage  here  is  that  there  is  no  need  for  the  merit  or  auxiliary  function  as  discussed  in  subsection  4.1  and  all 
function  values  and  subgradients  used  in  the  nonsmooth  maximization  process  are  computed  exactly  (i.e. , 
optimally).  Problem  decomposition  is  now  carried  out  for  the  n  dimensional  problem;  however,  decompo¬ 
sition  of  the  n  —  1  dimensional  recovery  problem  (and  all  lower  order  recovery  problems)  is  performed  only 
after  the  problem  is  formulated. 

Between  these  two  algorithms  are  a  host  of  different  relaxation  schemes  based  on  relaxing  n  —  m  sets 
of  constraints  to  an  m-dimensional  problem  (2  <  m  <  n),  but  these  all  have  the  same  difficulties  as  the  first 
algorithm  in  that  the  relaxed  problem  is  an  NP-hard  problem.  To  resolve  this  difficulty,  we  use  an  auxiliary 
or  merit  function  ^mn  as  described  in  subsection  4.2.  For  the  case  m  <  n  —  1,  the  recovery  procedure 
is  still  based  on  an  NP-hard  (n  —  m  +  l)-dimensional  assignment  problem.  The  decomposition  techniques 
discussed  in  subsection  4.7  are  based  on  identifying  the  assignment  problem  with  a  layered  graph  and  then 
finding  disjoint  components  of  this  graph.  In  general,  all  relaxed  problems  can  be  decomposed  prior  to 
any  nonsmooth  computations  because  their  structure  stays  fixed  throughout  the  algorithm.  All  recovery 
problems  cannot  be  decomposed  until  they  are  formulated,  as  their  structure  changes  as  the  solutions  to  the 
relaxed  problems  change. 

4.  Algorithm  Details  and  Refinements.  Many  aspects  of  the  algorithm  presented  in  Section  3  require 
further  explanation  on  how  various  problems  are  solved.  These  along  with  the  many  refinements  that  can 
significantly  increase  the  speed  of  the  relaxation  algorithm  are  explained  in  this  section. 

4.1  Maximization  of  the  Nonsmooth  Function  $mn(um+1 , . . . ,  un).  One  of  the  key  steps  in  the 
Lagrangian  relaxation  algorithm  in  subsection  3.3  is  the  solution  of  the  problem 

Maximize  ,  un)  |  uk  €  IRMfc+1  ;  k  =  m  +  1, . . . ,  n)  (4.1) 

where  =  0  for  ^  =  m  +  •  •  •  > n-  In  this  subsection  we  discuss  several  aspects  of  this  problem.  We 

first  show  that  this  is  a  problem  of  nonsmooth  optimization  and  then  briefly  explain  the  computation  of  the 
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required  subgradients.  Unfortunately,  the  evaluation  of  $mn  requires  the  evaluation  of  an  NP-hard  problem 
for  m  >  2.  Thus  for  real  time  needs,  we  address  this  difficulty  at  the  end  of  the  subsection.  The  actual 
nonsmooth  optimization  methods  are  presented  in  the  next  subsection. 

That  (4.1)  is  a  problem  of  nonsmooth  optimization  is  the  subject  of 

Theorem  4.1.  Let  um+1, . .  .  ,un  be  multiplier  vectors  associated  with  the  (m  +  l)5t  through  the  nth  set 
of  constraints  in  (3.1),  let  <£mn  be  as  defined  in  (3.2),  let  Vn(zn)  be  the  objective  function  value  of  the 
n-dimensional  assignment  problem  in  equation  (3.1),  let  zn  be  any  feasible  solution  of  (3.1),  and  let  zn 
be  an  optimal  solution  of  (3.1).  Then,  $mn(um+1 , . . . ,  un)  is  piecewise  affine,  concave  and  continuous  in 
{um+1, . . . ,  un }  and 

<W <  vn(zn)  <  Vn(zn),  (4.2 a) 

Furthermore, 

$mn(um+1 , . . . ,  un)  (4.26) 

for  all  uk  £  with  Uq  =  0  and  k  =  m, . . . ,  n. 

Thus  the  problem  of  maximizing  <£mn  is  one  of  nonsmooth  optimization. 

Most  of  the  algorithms  for  non-smooth  optimization  are  based  on  generalized  gradients  called  subgradi¬ 
ents,  given  by  the  following  definition. 

Definition  4.2.  At  u  =  (um+1 , . . . ,  un)  the  set  6$mn(u)  is  called  a  subdifferential  of  4>mn  and  is  defined  by 

=  {q  €  ]RMm+1+1  X  •  •  •  x  mM"  +  1  I  Smn(u>)  -  < 

qT{w -«)Vti)£  EM“+1+1  x  •  •  •  x  ]RM"+1}. 

A  vector  q  £  6$mn(u)  is  called  a  subgradient. 


If  zn  is  an  optimal  solution  of  (3.2)  computed  during  evaluation  of  $mn(«),  differentiating  $mn  with 
respect  to  u?n  yields  the  following  in-th  component  of  a  subgradient  g  of  <&mn(u) 

9o=  0 

Mi  Mk-i  Mk+ 1  Mn 

9ik  =Yi'"  X  X  12  2*v ~  1  for  ik  =  and  k  =  m+ 

*  1=0  tfe_i=0  ifc+1=0  *n=0 

If  2rn  is  the  unique  optimal  solution  of  (3.2),  6$mn(u)  =  { g },  and  $mn  is  differentiable  at  u.  If  the  op¬ 
timal  solution  is  not  unique,  then  there  are  finitely  many  such  solutions,  say  zn(l), ... ,  zn(K).  Given  the 
corresponding  subgradients,  g1, . . . ,  gK ,  the  sub  differential  6$(u)  is  the  convex  hull  of  {g1 , . . . ,  gK}  [24]. 
4.2  The  use  of  a  merit  or  auxiliary  function.  For  real-time  needs,  one  must  address  the  fact  that 
the  nonsmooth  optimization  problem  (4.1)  requires  the  solution  of  an  NP-hard  problem  for  m  >  2.  One 
approach  to  this  problem  is  to  use  the  following  merit  or  auxiliary  function  to  decide  whether  a  function 
value  has  increased  or  decreased  sufficiently  in  the  line  search  or  trust  region  methods: 

[  $mn(um+\...,un)  if  m  =  2 

^mn(«3) .  * . ,  um\ um+1 , . . . ,  un)  =  <  or  (3.2)  is  solved  optimally,  (4.5) 

[  $2n(^3) . . » j  um)  um+1, . . . ,  un )  otherwise. 

where  the  multipliers  u3,...,um  that  appear  in  lower  order  relaxations  used  to  construct  (subop timal) 
solutions  of  the  m-dimensional  relaxed  problem  (3.2)  have  been  explicitly  included.  Note  that  ^  is  well- 
defined  since  (3.4)  can  always  be  solved  optimally  if  m  =  2.  For  sufficiently  small  problems  (3.2)  or  (3.4),  one 
can  more  efficiently  solve  the  NP-hard  problem  by  branch  and  bound.  This  is  the  reason  for  the  inclusion 
of  the  first  case;  otherwise,  the  relaxed  function  $2 n  to  guide  the  nonsmooth  optimization  phase.  That  the 
merit  function  provides  a  lower  bound  for  the  optimal  solution  follows  directly  from  Theorem  4.1  and 
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Theorem  4.3.  Given  the  definition  of  ^ mn  in  (4.5) 


vmn{u3 


,um:um+\ 


u")  <  $mn(wm+1, . . . ,  tx”)  for  all  multipliers  u3, . . . ,  um,  um+1, . 


un.  (4.6) 


The  actual  function  value  used  in  the  optimization  phase  is  lmn;  however,  the  subgradients  are  computed 
as  in  (4.4),  but  with  the  solution  in  being  a  suboptimal  solution  constructed  from  a  relaxation  procedure 
applied  to  the  m-dimensional  problem.  Again,  the  use  of  these  lower  order  relaxed  problems  is  the  reason 
for  the  inclusion  of  the  multipliers  u3, . . . ,  um . 

To  explain  how  the  merit  function  is  used,  suppose  we  have  a  current  multiplier  set  (u™^1, . . .  ,u£ld) 
and  we  wish  to  update  to  a  new  multiplier  set  (u™*1  >  ■  •  • »  u£ew)  via  Wet1 1  •  •  •  >  Wnew)  =  Wd  1  >  •  •  •  >  w2id)  + 
(Aitm+1, . . . ,  Aun).  Then  we  compute  Vmn  (u3ld, ...  ,u^d;  u"ld)  where  (u3ld,  *  ■  •  Wd)  is  obtained 

during  the  relaxation  process  used  to  compute  a  high  quality  solution  to  the  relaxed  m-dimensional  assign¬ 
ment  problem  (3.2)  at  (um+1, . .  .,un)  =  (ti^J1,  .  -  •  > «Sld)-  The  decision  to  accept  (it™*1,  •  •  .,«Sew)  is  then 
based  on  mn  (u3  d , . . . ,  u™+1 , . . . ,  u*ld)  <  (<4W *  •  •  •  >  CewJ  Cet1 ,  ■  •  ■ ,  <ew)  or  some  other  stopping 

criteria  commonly  used  in  line  searches.  Again,  (u3ew, . . . ,  u™ew)  represents  the  multiplier  set  used  in  the  lower 
level  relaxation  procedure  to  construct  a  high  quality  feasible  solution  to  the  m-dimensional  relaxed  problem 
(3.2)  at  (um+1, . . . ,  un)  =  (iCt1  >  •  •  • » «new )■  P°int  is  that  each  time  one  changes  (um+1 , . . . ,  un)  and 
uses  the  merit  function  ^mn(u3, . . . ,  um]  um+ 1 , . . . ,  un)  for  comparison  purposes,  one  must  generally  change 
the  lower  level  multipliers  (u3, . . . ,  um). 

An  illustration  of  this  merit  function  for  m  =  n  -  1  is  given  in  the  work  of  Poore  and  Rijavec  [50]. 

4.3.  Nonsmooth  Optimization  Methods.  By  Theorem  4.1  the  function  4>mn(u)  is  a  continuous, 
piecewise  affine,  and  concave,  so  that  the  negative  of  $mn(u)  is  convex.  Thus  the  problem  of  maximizing 
$mn(u)  is  one  of  nonsmooth  optimization.  There  is  a  large  amount  of  literature  on  such  problems  [29,  30, 
33,  62,  63,  65].  Suffice  it  to  say  that  we  have  tried  a  variety  of  methods  including  subgradient  methods 
[63],  bundle  methods  [29,  30,  33],  and  the  recent  bundle  trust  method  of  Schramm  and  Zowe  [62].  We  have 
determined  that  for  a  fixed  number  of  nonsmooth  iterations,  say,  ten,  the  bundle-trust  method  provides  good 
quality  solutions  with  the  fewest  number  of  function  and  subgradient  evaluations  of  all  the  methods,  and  is 
therefore  our  currently  recommended  approach. 

4.4.  The  Two  Dimensional  Assignment  Problem.  The  forward/reverse  auction  algorithm  of  Bert- 
sekas,  Castanon,  and  Tsaknakis  [9]  is  used  to  solve  the  many  relaxed  two  dimensional  problems  that  occur 
in  the  course  of  execution. 

4.5.  Initial  Multipliers  and  Hot  Starts.  The  effective  use  of  “hot  starts”  is  fundamental  for  real-time 
applications.  A  good  initial  set  of  multipliers  can  significantly  reduce  the  number  of  nonsmooth  iterations 
(and  hence  the  number  of  4>mn  evaluations)  required  for  a  high  quality  recovered  solution.  A  presentation 
of  these  techniques  can  be  found  in  the  thesis  of  Robertson  [60]. 

4.6  Local  Search  Methods.  Given  a  feasible  solution  of  the  multidimensional  assignment  problem,  one 
can  consider  local  search  procedures  to  improve  this  result  [38].  A  discussion  of  these  methods  is  presented 
in  the  thesis  of  Robertson  [60]. 

4.7.  Problem  Decomposition.  The  algorithm  described  thus  far  is  all  based  on  relaxation.  Due  to  the 
sparsity  of  the  problems,  one  can  frequently  decompose  the  problem  into  a  collection  of  disjoint  components 
each  of  which  can  be  solved  independently.  Due  to  the  setup  costs  of  Lagrangian  relaxation,  a  branch  and 
bound  procedure  is  generally  more  efficient  for  small  components,  say  four  or  five  feasible  arcs.  Otherwise, 
we  use  the  relaxation  procedures  described  above. 
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Perhaps  the  easiest  way  to  view  our  decomposition  method  is  to  view  the  reports  or  measurements  as 
a  layered  graph.  A  vertex  is  associated  with  each  observation  point,  and  an  edge  is  allowed  to  connect  two 
vertices  only  if  the  two  observations  belong  to  at  least  one  feasible  track  of  observations.  Given  this  graph, 
the  decomposition  problem  can  then  be  posed  as  that  of  identifying  the  connected  subcomponents  of  a  graph 
which  can  be  accomplished  by  constructing  a  spanning  forest  via  a  depth  first  search  algorithm  [1]. 

The  orginal  relaxation  problem  is  decomposed  first.  All  relaxed  assignment  problems  can  be  decomposed 
a  priori  and  all  recovery  problems  can  be  decomposed  only  after  they  are  formulated.  Hence,  in  the  n-to- 
(n  -  1)  case,  we  have  n  -  2  relaxed  problems  that  can  all  be  decomposed  initially,  and  the  recovery  problems 
are  not  decomposed  (since  they  are  all  two  dimensional).  In  the  n-to-2  case,  we  have  only  one  relaxed  problem 
that  can  be  decomposed  initially.  This  case  yields  n  -  3  recovery  problems,  which  can  be  decomposed  only 
after  they  are  formulated. 

4.8  Use  of  Invariance  and  Cost  Shifting.  A  heuristic  that  can  add  significant  speed  to  the  overall 
relaxation  procedure  is  to  use  the  following  invariance  principle  to  shift  the  most  negative  cost  to  zero .  This 
work  is  discussed  in  the  work  of  Rijavec  and  Poore  [50]. 

5.  Numerical  Performance  of  the  Algorithm.  The  Lagrangian  relaxation  algorithm  in  subsection 
3.3  with  m  =  2  is  sufficiently  general  to  encompass  a  wide  range  of  dense  and  sparse  multidimensional 
assignment  problems.  We  should  note  a  few  features  of  our  implementation  that  have  an  impact  on  the 
performance  numbers  that  follow.  First,  the  problems  considered  do  not  decompose.  Second,  we  store 
only  free  variables  (as  opposed  to  a  multidimensional  matrix).  The  algorithm  incorporates  the  bundle-trust 
method  [62]  to  solve  the  nonsmooth  optimization  problem  (4.1)  and  the  forward/reverse  asymmetric  auction 
algorithm  [8,  9]  to  solve  the  relaxed  two-dimensional  assignment  problem  (3.4).  The  program  currently 
executes  the  same  number  of  nonsmooth  optimization  steps  for  each  assignment  problem  for  which  n  >  2. 
When  the  solver  reaches  a  3-dimensional  recovery  problem,  a  2-dimensional  recovery  procedure  is  executed 
for  each  nonsmooth  iteration. 

As  we  now  demonstrate  for  a  particular  class  of  problems,  the  execution  time  for  the  algorithm  in 
subsection  3.3  with  m  =  2  is  linear  in  the  number  of  free  or  feasible  variables  and  in  the  number  of  nonsmooth 
iterations.  This  is  due  to  the  fact  that  the  dominant  computational  part  of  the  algorithm  is  the  evaluation 
of  the  relaxed  two-dimensional  cost  coefficients  as  specified  in  (3.3).  The  problem  class  that  we  have  chosen 
is  4-dimensional  with  Mk  =  25  on  average  for  k  =  1, . . . ,  4.  The  cost  coefficients  are  uniformly  random  on 
the  interval  [-100,-1].  The  times  are  averages  over  100  randomly  generated  problems,  with  the  algorithm 
running  on  an  IBM  RISC  6000/550.  Our  first  test  was  to  determine  execution  time  as  a  function  of  the 
number  of  free  variables  for  a  fixed  number  of  nonsmooth  iterations.  Our  second  test  was  to  determine 
execution  time  as  a  function  of  the  number  of  nonsmooth  iterations,  with  a  fixed  number  of  free  variables. 
Results  are  shown  in  Figures  5.1  and  5.2. 

In  these  graphs,  the  lines  correspond  to  a  linear  least  squares  fit  of  the  data  provided.  In  Figure  5.1,  the 
upper  line  represents  the  overall  execution  time  and  the  lower  line  represents  the  time  spent  in  evaluation  of 
the  relaxed  cost  coefficients.  Figure  5.1  shows  that  the  growth  in  execution  time  is  a  linear  function  of  the 
number  of  free  variables,  and  the  evaluation  of  the  relaxed  cost  coefficients  (3.3)  clearly  dominates  overall  ex¬ 
ecution.  Also,  the  fact  that  the  gap  grows  slightly  between  the  two  lines  as  the  number  of  arcs  increases  shows 
that  the  other  algorithm  components  grow  slightly  in  execution  time.  Figure  5.2  shows  that  the  growth  in  ex¬ 
ecution  is  a  linear  function  of  the  number  of  nonsmooth  iterations.  We  can  therefore  conclude  that  the  overall 
performance  time  of  algorithm  in  subsection  3.3  with  m  =  2  is  a  linear  function  of  the  number  of  free  variables. 
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•Total  Time 
•Afg  Min  Time 


Number  of  Arcs 


Figure  5.1  -  Number  of  Arcs  vs  Speed 


Figure  5.2  -  Number  of  Nonsmooth  Iterations  vs  Speed 


Next  we  consider  the  comparison  of  two  different  suboptimal  methods  (relaxation  and  randomized 
greedy)  to  the  optimal  solution  provided  by  branch- and-bound,  run  for  50  randomly  generated  problems. 
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The  randomized  greedy  was  the  best  of  serveral  other  heuristic  methods  including  greedy  and  max-regret, 
and  thus  we  display  on  the  results  for  this  randomized  greedy.  The  problems  were  all  four  dimensional, 
with  average  size  M*  =  7,  and  costs  uniformly  generated  on  the  interval  [1, 100].  The  relaxation  method 
was  executed  to  convergence,  which  averaged  200  nonsmooth  iterations.  All  objective  function  values  are 
normalized  to  the  optimal  (100)  solution.  The  graph  shows  the  superior  quality  of  the  recovered  feasible 
solutions,  and  also  demonstrates  the  small  approximate  duality  gaps  provided  by  relaxation. 


6.  Tracking  Software.  Since  the  beginning  of  1994,  we  have  completely  redesigned  and  rewritten  our 
1993  multitarget  tracker  to  include  multisensor  processing  with  new  and  improved  data  structures.  The 
current  program,  which  is  about  twenty  thousand  lines  of  C  code,  is  modular  and  is  designed  for  adaptation 
to  different  sensor  fusion  scenarios:  moving  vs.  stationary  platforms,  collocated  vs.  distributed  sensors, 
passive  and/or  active  sensors,  and  handles  asynchronous,  out  of  sequence  and  delayed  measurements/data. 
New  maneuvering  target  algorithms  now  utilize  the  unique  capability  for  fast  re-association  over  multiscan 
window  to  switch  between  multiple  models  for  the  target  dynamics.  This  software  is,  however,  not  available 
at  the  time  this  report  was  written. 

7.  A  Two  Radar  System  with  Rome  Labs  Data.  This  section  describes  a  tracking  problem  based  on 
six  and  a  half  minutes  of  real  data  observed  by  radars  at  Dansville  and  Remsen  in  upstate  New  York.  The 
Remsen  radar  is  approximately  20  NM  Northeast  of  Rome  Labs  while  the  Dansville  radar  resides  about  100 
NM  West  Southwest  of  Rome  Labs.  Both  radars  are  L-band  radars.  They  consist  of  a  search  radar  and  a 
Beacon  Interrogator  with  a  common  digitizer  set.  The  antenna  of  the  radar  and  the  beacon  are  collocated 
and  move  together.  The  JSS  site  characteristics  of  interest  are: 
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Antenna  speed 

Frequency 

Prf 

Pulse  width 
Vertical  beamwidth 
Horizontal  beamwidth 
Range 


DANSVILLE,  NY 
6  rpm 

1280-1350  MHz 
350  Hz 
1.8  microsec 
3.75° 

1.2° 

200  NM 


REMSEN,  NY 

5  rpm 

1250  -  1350  MHz 
340  Hz 

6  microsec 
5.4° 

1.3° 

200  NM 


The  data  is  received  by  specialized  hardware  called  Data  Interface  Unit  (DIU)  where  the  asynchronous 
data  is  formatted  into  scans  from  each  radar  and  sent  on  to  files  on  a  Micro  VAX  II  or  used  as  input  to 
fusion  algorithms.  The  detection  times  are  generated  by  the  combining  of  Universal  Time  code  (UTC)  with 
the  delta  time  supplied  by  each  input  detection  message.  The  two  types  of  sensors  collocated  at  the  radar 
site  detect  two  types  of  targets:  radar  and  beacon.  The  radar  system  picks  up  the  reflected  signal,  or  “skin 
paint,”  from  a  non-cooperative  target.  Note  that  the  signal  can  also  be  reflected  by  weather,  atmospheric 
conditions,  terrain,  etc.,  causing  false  targets  or  clutter.  Various  canceling  mechanisms  can  remove  much  of 
this  clutter,  but  the  remainder,  or  “clutter  residue,”  appears  to  the  system  as  radar  targets.  The  Beacon  IFF 
(Identification  Friend  of  Foe)  system  relies  on  a  transponder  onboard  the  cooperative  aircraft  to  respond  to 
“who  are  you?”  and  “how  high  are  you?”  queries  from  the  site.  The  Beacon  system  is  thus  of  limited  utility 
in  a  threat  scenario,  except  to  identify  targets  that  are  probably  not  a  threat.  It  is  also  useful  for  providing 
“ground  truth”  for  system  testing  when  a  cooperative  target  is  used  for  a  simulated  threat. 

The  beacon  responds  with  an  IFF  response  (squawk)  and  the  altitude  of  the  target.  This  response  is 
a  four  digit  number  that  is  set  manually  by  the  target  aircraft  denoting  the  identity  of  the  aircraft.  This 
number  may  change  during  flight  under  direction  of  air  traffic  control.  If  either  of  the  last  two  digits  are  non 
zero,  the  code  is  said  to  be  discrete  and  unique  to  one  aircraft  in  the  sky.  Usually  all  aircraft  that  are  flying 
visual  flight  rules  (VFR)  are  assigned  the  same  number,  1200.  Some  aircraft  do  not  have  transponders.  The 
altitude  of  the  aircraft  is  sent  with  the  beacon  response.  The  altitude  is  a  barometric  pressure  reading.  A 
barometric  altitude  correction  is  needed  for  altitudes  below  18,000  feet  to  correct  it  to  true  feet  above  Mean 
Sea  Level.  An  apparent  discontinuity  in  the  altitude  could  occur  as  the  target  climbs  or  descends  through 
the  “transition  altitude”  of  18,000  feet.  In  summary,  there  are  three  types  of  target  detections  reported: 

Search  Radar  Detection:  This  is  a  detection  by  the  radar  system  only.  The  target  could  be  oriented  such 
that  the  beacon  signal  did  not  activate  the  target’s  beacon  or  the  target’s  beacon  system  could  be  turned 
off,  but  the  target  returned  a  radar  signal  (skin  return).  This  detection  could  also  be  a  false  alarm. 

Beacon  Reinforced  Target:  This  detection  report  results  when  both  the  beacon  receives  a  transponder 
return  and  the  radar  receives  a  skin  return  at  the  same  position. 

Beacon  Only  Target:  A  signal  received  from  the  target’s  beacon  which  is  not  correlated  with  a  radar 
signal.  This  can  occur  when  the  target  aircraft’s  radar  cross  section  is  such  that  the  radar  skin  return  is 
below  the  detection  threshold,  but  the  one  way  energy  of  the  beacon  signal  is  strong  enough  to  trigger  the 
target  aircraft’s  transponder. 

More  than  sixty  aircraft  are  observed  in  the  period  covered  by  the  data.  There  is  also  some  clutter, 
both  random  and  stationary,  probably  due  to  ground  features.  Almost  all  the  observations  in  the  problem 
come  from  the  overlap  region(i.e.,  the  region  observed  by  both  radars).  Several  aircraft,  however,  are  shown 
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outside  of  that  area,  especially  north  of  Remsen.  Due  to  the  terrain  features,  and  possibly  other  causes,  a 
number  of  targets  within  the  overlap  region  are  observed  by  one  radar  only  through  all  or  part  of  their  track 
life.  For  some  targets,  even  the  coverage  by  one  sensor  is  intermittent. 

Rome  Data  -  Flat  Earth  (All) 


Figure  7.1. 

At  the  ranges  present  in  the  problem,  the  curvature  of  the  Earth  becomes  an  issue  and  the  “flat  earth 
assumption”  leads  to  significant  misalignment  problems.  Figure  7.1  shows  the  observations  in  the  tracking 
plane  generated  by  ignoring  the  Earth’s  curvature.  It  can  be  seen  that  the  observations  from  the  two  sensors 
form  parallel  lines  for  each  target  observed  by  both  sensors. 
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Figure  7.2  shows  that  the  curvature  of  the  Earth  must  be  taken  into  account.  Mulholland  and  Stout  [36] 
describe  a  stereographic  projection  algorithm  used  in  the  National  Aerospace  System.  The  tradeoffs  between 
the  different  projections  are  described  in  some  detail  by  Goldenberg  and  Wolf,  [25]  and  the  stereographic 
projection  is  determined  to  be  superior  to  other  alternatives.  The  stereographic  projection  was  thus  chosen 
to  transform  the  tracking  problem  into  a  Cartesian  coordinate  system. 

Aligned  Rome  Data  (All) 


Figure  7.2. 

The  tracking  coordinate  system  was  chosen  to  be  3D  Cartesian  coordinate  system,  with  the  origin 
midpoint  between  the  radar  sites,  and  the  height  as  average  of  the  height  of  the  radar  sites.  The  positive  x 
axis  pointed  east,  while  the  positive  y  axis  pointed  north. 
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Even  after  applying  the  stereographic  projection,  a  slight  misalignment  remained  in  the  data.  It  was 
found,  however,  that  applying  a  0.5  kilometer  correction  to  the  range  of  the  Remsen  radar  and  0.1°  correction 
to  the  azimuth  of  the  Dansville  radar  removed  almost  all  the  misalignments.  The  aligned  observations  are 
shown  in  Figure  7.3.  When  more  data  becomes  available,  a  more  thorough  analysis  of  the  alignment  problems 
can  be  made  using  system  identification  algorithms. 

Computing  a  stereographic  projection  of  an  observation  with  known  height  is  straightforward.  If,  how¬ 
ever,  only  range  and  azimuth  are  available,  the  value  of  the  height  parameter  must  be  assumed.  This  was 
taken  to  be  18,000  feet.  The  additional  inaccuracy  in  the  data  resulting  from  the  height  assumption  must 
be  addressed  in  the  tracking  algorithms. 

7.2  Numerical  Solution.  This  section  will  give  a  brief  discussion  of  the  algorithms  used  and  experience 
gained  in  solving  the  two  radar  tracking  problem  described  in  the  preceding  section.  In  a  general  multisensor 
tracking  problem,  the  target  space  would  likely  be  partitioned  according  to  the  sensor  coverage,  and  each 
part  tracked  separately,  with  additional  logic  to  handle  tracks  crossing  from  one  region  into  another.  Since 
almost  all  the  observations  in  this  particular  problem  lay  in  the  sensor  overlap  region,  such  partitioning  was 
not  necessary. 

The  data  stream  from  each  radar  was  partitioned  into  scans.  Since  both  sensors  made  sweeps,  as  opposed 
to  instantaneous  snapshots,  of  space,  some  observations  were  perceived  by  the  tracker  as  out  of  time  order. 
This  was  somewhat  exacerbated  by  the  fact  that  the  Dansville  radar  has  a  faster  scan  rate  than  the  Remsen 
radar.  As  a  result,  two  scans  from  Dansville  radar  would  occasionally  be  sequenced  one  after  another  for 
processing.  This  meant  up  to  three  observations  could  be  out  of  time  order.  Out  of  sequence  measurement 
problem  would  become  more  severe  if  more  sensors  were  used  or  the  difference  in  scan  rates  was  larger. 

Handling  out  of  sequence  measurements  imposes  a  significant  programming  overhead  on  the  tracking 
software,  even  though  it  does  not  involve  a  measurable  run-time  penalty.  The  complications  arise  in  gating, 
estimation,  maneuver  detection  and  output.  Our  current  tracker  is  designed  to  handle  sensors  with  widely 
varying  scan  rates,  resulting  in  almost  arbitrary  number  of  measurements  that  arrive  out  of  sequence. 

To  compute  the  target  tracks,  an  iterated  extended  Kalman  filter  with  multiple  models  was  used  for  the 
different  target  dynamics.  Figure  7.3  presents  the  tracks  that  were  constructed  by  the  tracker  as  a  solution 
of  the  tracking  problem  shown  in  Figure  7.2.  No  additional  smoothing  was  done  for  the  tracks  being  output. 

Height  estimation  for  the  non-transponder  targets  poses  an  additional  challenge.  Numerical  experience 
using  the  data  generated  by  a  simulator  indicates  that  a  system  identification  algorithm  can  estimate  the 
target  height  very  accurately  if  the  target  is  observed  by  two  radars.  The  accuracy  is  typically  within  a  few 
hundred  meters,  even  though  a  simple  triangulation  could  possibly  lead  to  accuracy  of  no  more  than  3000 
meters  given  the  range  and  azimuth  measurement  errors. 

In  handling  the  real-world  data,  the  stereographic  projection  algorithm  imposes  additional  inaccuracy, 
since  a  wrong  assumed  height  of  the  target  will  result  in  an  incorrect  projection.  If  the  transponder-equipped 
targets  in  the  problem  (whose  height  is  known)  are  treated  as  if  their  height  is  not  known,  the  preliminary 
numerical  experience  shows  that  the  height  can  be  estimated,  projection  errors  notwithstanding.  Even  in 
the  cases  where  the  targets  were  observed  by  mostly  one  radar,  with  only  an  occasional  observation  by  the 
second  radar,  the  system  identification  algorithm  gave  a  height  estimate  accurate  to  within  1200  meters. 
Due  to  the  limited  data  sample,  however,  the  final  analysis  of  the  height  estimation  algorithms  will  have  to 
wait  until  more  data  becomes  available. 
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Rome  Data  -  Full  Problem  Solution 


Figure  7.3. 
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3.  Data  association  problems  posed  as  multidimensional  assignment  problems:  problem  formulation,  SPIE, 
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